C.............................. Get_Data.FOR ...............................
C...... This file contains routines to supply data. These routines may be replaced by the user.
C:::::::::::::::::::::::::::::::: GET_EUV ::::::::::::::::::::::
C....... This subroutine just supplies a sample SEE EUV spectrum
C....... P. Richards June 2015
      SUBROUTINE GET_EUV(FDIM,  !.. IN: EUV Flux dimension
     >                   NEUV,  !.. OUT: # of EUV fluxes
     >                EUVFLUX,  !.. OUT: EUV Fluxes (Watts)
     >                   WLNM)  !.. OUT: EUV wavelengths (nm)
      IMPLICIT NONE
      INTEGER I                  !.. Looping variable
      INTEGER FDIM,NEUV          !.. See I/O parameters
      INTEGER NFLUXDAT           !.. # of solar EUV fluxes in DATA block
      PARAMETER (NFLUXDAT=100)
      REAL EUVDAT(NFLUXDAT),EUVFLUX(FDIM),WLNM(FDIM)

c      !.. SEE solar irradiances (x1.0E7) for June 2007 at 1 nm intervals
c      DATA EUVDAT/18,312,188,96,236,204,308,197,148,147,105,89,58,91,58,
c     > 391,826,2107,1351,798,694,744,443,259,323,501,201,914,862,765,
c     > 2821,615,328,317,390,285,381,125,69,48,67,48,47,98,56,55,151,
c     > 65,98,119,134,51,56,64,39,223,52,44,337,66,130,104,279,231,34,
c     > 30,34,27,61,42,96,38,33,24,35,49,128,110,159,106,83,88,105,
c     > 235,144,172,193,227,263,296,338,311,118,123,93,72,55,734,115,156/

      !.. New SEE solar irradiances (x1.0E7) for June 28, 2007 at 1 nm
      DATA EUVDAT/17,319,209,136,238,218,238,313,293,264,198,147,77,85,
     > 169,132,386,2152,1160,805,472,407,608,368,428,677,166,584,757,
     > 707,3046,620,314,370,418,318,436,126,64,42,60,46,42,92,51,54,
     > 150,64,96,120,134,53,58,64,39,218,52,45,337,64,133,106,273,228,
     > 34,31,33,27,61,42,96,37,33,24,35,48,129,110,158,106,82,87,102,
     > 244,141,174,193,230,266,301,351,320,118,124,93,72,57,756,118,161/
     
      NEUV=NFLUXDAT   !.. # solar EUV fluxes to return

      !.. Fill up the EUV flux (Watts) and wavelength (nm) arrays for PE2S
      WLNM(1)=0.5
      DO I=1,NFLUXDAT
        EUVFLUX(I)=EUVDAT(I)*1.0E-7
        IF(I.GT.1) WLNM(I)=WLNM(I-1)+1.0
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::::::::::::: GET_EDEN :::::::::::::::::::::::::::::::
C...... This routine returns electron density from a profile based on NmF2 and hmF2
C...... P.Richards, June 2015
      SUBROUTINE GET_EDEN(ALT,HMF2,NMF2,ZNE)
      IMPLICIT NONE
      REAL ALT,HMF2,NMF2,ZNE
      REAL HLBY,NLBY                !.. Lower boundary height and density
      REAL F2CHAP,SHF2,F2EXP,SHEXP  !.. F2 region parameters
      DATA HLBY, SHF2, SHEXP/100., 70.0, 120.0/

      !.. The lower boundary density variation is based on NmF2 (solar cycle)
      NLBY=SQRT(NMF2/2.0E5)*7E4
      !.. Linear variation between 100 km and hmF2
      ZNE=((HMF2-ALT)*NLBY+(ALT-HLBY)*NMF2)/(HMF2-HLBY)
      !.. Chapman variation near the F2 peak
      F2CHAP=NMF2*(1-(HMF2-ALT)*(HMF2-ALT)/4/SHF2/SHF2)
      !.. Below hmF2, density is the maximum of current Ne and F2 Chapman
      IF(ALT.LT.HMF2) THEN
          ZNE=MAX(ZNE,F2CHAP)
      ELSE
         ZNE=F2CHAP
      ENDIF
      !.. Above hmF2, density is the maximum of current value and an
      !.. F2 exponential variation
      F2EXP=0.0
      IF(ALT.GT.HMF2) THEN
        F2EXP=NMF2*EXP((HMF2-ALT)/SHEXP)       
        ZNE=(ZNE+F2EXP)/2.0
        ZNE=MAX(ZNE,F2EXP)
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::: BRACE_TE :::::::::::::::::::::::::::::::
C...... Empirical Te model from Brace and Theis GRL p275, 1978
      SUBROUTINE BRACE_TE(H,NI,TN,TB)
      IMPLICIT NONE
      REAL H,NI,TN,TB,TCON
      TB=1051+(17.07*H-2746)*EXP(-5.122E-4*H+6.094E-6*NI-3.353E-8*H*NI)
      IF(TB.LT.TN) TB=TN
      TCON=(TB+TN)/2.0
      !.. Restrain Te at high altitudes where BRACE does not apply
      IF(TB.GT.2500.0) TB=(TB+2500.0)/2.0
      RETURN
      END
C::::::::::::::::::::::::::::::::::::  :::::::::::::::::::::::::::::
C..... This subroutine reads in the F10.7 and Kp indices from the master file,
C..... converts Kp to Ap and returns AP and F10.7. 
C..... Written by P. Richards March 2009. Updated 2015
      SUBROUTINE READ_F10KP(INDATE,      !.. Date (YYYYddd)
     >                        F107,      !.. OUT: F10.7 
     >                       F107A,      !.. OUT: F10.7A 
     >                       APDAY)      !.. OUT: Daily Ap
      IMPLICIT NONE
      INTEGER I,J         !.. loop control variables
      INTEGER INDATE
      INTEGER YYYY,MM,DD,DDD,YYYYDDD,KP(8),KPTOAP(0:90)
      REAL KPDAY,KPAVE,F107,F107A,APDAY

      !.. Conversion factors from Kp = 1..9. AP = KPTOAP(KP)
      !.. Kp= 0    1    2   3    4    5    6     7     8    9
      !.. AP= 0,   3,   7,  15,  27,  48,  80,  140,  240, 400
      !.. 90 Ap values for easy conversion to Ap
      DATA KPTOAP/3*0,4*2,3*3,3*4,4*5,3*6,3*7,4*9,3*12,3*15,4*18,3*22
     >   ,3*27,4*32,3*39,3*49,4*56,3*67,3*80,4*94,3*111,3*132,4*154
     >   ,3*179,3*207,4*236,3*300,400/

      !.. read headers off the F10.7 and Kp input file
 10   READ(20,*,ERR=10,END=99) YYYY,MM,DD,(KP(J),J=1,8),F107,F107A

      !.. Loop to find the specified date
      DO I=1,9999999
        KPAVE=REAL((KP(1)+KP(2)+KP(3)+KP(4)+KP(5)+KP(6)+KP(7)+KP(8))/8)
        APDAY=KPTOAP(NINT(KPAVE))
        KPDAY=KPAVE/10.0
        CALL CONVT_DATE(MM,DD,DDD,YYYY)  !.. Convert date to DDD
        YYYYDDD=YYYY*1000+DDD 
        IF(INDATE.LT.YYYYDDD) THEN
          WRITE(6,*) 
     >     ' ** Start date is before first date in F10Kp file **'
          STOP
        ENDIF
        IF(INDATE.EQ.YYYYDDD) RETURN
        READ(20,*,ERR=99,END=99) YYYY,MM,DD,(KP(J),J=1,8),F107,F107A
      ENDDO
      RETURN
 99   CONTINUE
      WRITE(6,*) ' ** Not enough days in the F10Kp file  **'
      STOP
      END
C::::::::::::::::::::::::: SUBROUTINE CONVT_DATE :::::::::::::::::::::::::::
C---- Scot A Braze,  8/7/95
C---- This subroutine converts the date from YYYY MM DD to YYYY DDD. 
      SUBROUTINE CONVT_DATE(MM,DD,DDD,YYYY)
      IMPLICIT NONE
      INTEGER DDD,MM,DD,NCONV(12),LYCONV(12),YYYY

      DATA NCONV/0,31,59,90,120,151,181,212,243,273,304,334/
      DATA LYCONV/0,31,60,91,121,152,182,213,244,274,305,335/

      IF(MOD(YYYY,4).EQ.0) THEN
          DDD=LYCONV(MM)+DD
      ELSE
          DDD=NCONV(MM)+DD
      ENDIF
      END
C::::::::::::::::::::::::::::::::::::: GETSZA ::::::::::::::::::::::::::::::::
C...... Routine to calculate solar zenith angle and local time
      SUBROUTINE GETSZA(YYYYDDD, !.. IN: Date
     >                       UT, !.. IN: UT in hours
     >                    GLATD, !.. IN: Latitude (degrees)
     >                    GLOND, !.. IN: Longitude (degrees)
     >                      SAT, !.. OUT: Solar apparent time (hours)
     >                      SZA) !.. OUT: solar zenith angle (degrees)
      IMPLICIT NONE
      INTEGER YYYYDDD
      REAL UT,GLATD,GLOND,SAT,SZA
      REAL DECL,ETRAN,SZA_ARG,HH,GLAT
      
      CALL SOLDEC(YYYYDDD,UT,DECL,ETRAN)  !.. solar declination, ephemeris transit

      !...... Evaluate Solar apparent time
      SAT=(UT*3600-ETRAN+43200.0)/3600+GLOND/15.0
      IF(SAT.LT.0) SAT=SAT+24
      IF(SAT.GT.24) SAT=SAT-24

      !...... Evaluate Solar zenith angle
      HH=(SAT-12.)*15.0*3.1415926/180.       !.. hour angle
      GLAT=GLATD/57.29578
      !.. Evaluate argument for ACOS and test because roundoff error 
      !.. may cause invalid ACOS argument      
      SZA_ARG=COS(GLAT)*COS(DECL)*COS(HH)+SIN(GLAT)*SIN(DECL)
      IF(SZA_ARG.GT.0.999) SZA_ARG=0.999
      IF(SZA_ARG.LT.-0.999) SZA_ARG=-0.999
      SZA=ACOS(SZA_ARG)*57.2958               !.. solar zenith angle
      RETURN
      END
C::::::::::::::::::::::::::: SOLDEC :::::::::::::::::::::::::
C...... This routine calculates the solar declination (DELTA radians)
C...... and ephemeris transit time (ETRAN in secs). ETRAN is the time
C...... the sun is overhead in Greenwich
C...... REFERENCE - page 484 "Explanatory Supplement to the Astronomical
C...... Almanac" Kenneth Seidelmann, University Science Books, 20 Edgehill 
C...... Road, Mill Valley, CA 94941, 1992
C...... IDAY is yyyyddd (eg 1988015 = feb 15), UT is the universal time
C...... in hours (0-24) 
       SUBROUTINE SOLDEC(IDAY,UT,DELTA,ETRAN)
      !.... Find day of year (0-366), and year (1996), then Julian day
       INTEGER JD,DAYNUM
       REAL T,YEAR,UT,L,G,LAMBDA,EPSIL,E,GHA,DELTA,SD,ETRAN,DTOR
       DOUBLE PRECISION DJD,DUT
       DATA DTOR/57.29578/   ! degrees to radians
      !..... Recover year and date and make sure UT is less than 24
       DAYNUM=MOD(IDAY,1000)
       YEAR=(IDAY-DAYNUM)/1000
       JD=INT(365.25*(YEAR-1900)+DAYNUM)+2415020      ! Julian day
       DJD=JD                                ! extra precision needed         
       DUT=UT
       T=(DJD+DUT/24.0-2451545.0)/36525.0         ! # of centuries
       L=AMOD(280.460+36000.770*T,360.0)              ! aberration
       G=AMOD(357.528+35999.050*T,360.0)              ! mean anomaly
      !...... LAMBDA= ecliptic longitude. DELTA=obliquity of the ecliptic.
       LAMBDA=AMOD(L+1.915*SIN(G/DTOR)+0.020*SIN(2.0*G/DTOR),360.0)
       EPSIL=23.4393-0.01300*T
      !...... Equation of time. Time difference between noon and overhead sun
       E=-1.915*SIN(G/DTOR)-0.020*SIN(2.0*G/DTOR)
     >   +2.466*SIN(2*LAMBDA/DTOR)-0.053*SIN(4*LAMBDA/DTOR)
       GHA=15*UT-180+E                    ! Greenwich hour angle
       DELTA=ASIN(SIN(EPSIL/DTOR)*SIN(LAMBDA/DTOR))  ! solar declination
       SD=0.267/(1-0.017*COS(G/DTOR))    ! 
       ETRAN=(12-E/15)*3600 ! Ephemeris transit time in secs for FLIP
       RETURN
       END

